Protocol 1 - Geochemical data
The following example applies protocol 1 to confirm the workshops’ chemical reference groups.
Protocol 1 consist in:
- Select geochemical compositional data (CoDa);
- Apply transformation;
- Perform robust Principal Components Analysis (robPCA), implicitly using Euclidean distance;
- Perform PERMANOVA & PERMDISP tests;
Last, search for outliers and re-do protocol excluding outliers.
NOTE: The initial procedures must be ran at least once before any protocol can be applied.
Ordination procedure
prot1 <- apply_ordination(cleanAmphorae[!isShipwreck,], # no shipwrecks
"1", # select protocol 1
coda_override = chemVars16,
coda_transformation = "ILR")
#> sROC 0.1-2 loaded
#> [1] "78.24% of variance explained in 2D"
#> [1] "86.54% of variance explained in 3D"
#> [1] "Protocol 1 ended."The outcome is an ordination object. In this case, it is the output of pcaCoDa function in robCompositions package, in addition to several extra information, such as the transformed data (‘transformed_data’), the distance matrix (‘dist_matrix’), and the ready-to-plot texts indicating the fitness of the 2D/3D projections respect the distance matrix (‘sub2d’, ‘sub3d’). The later are printed in the console once the object is created.
class(prot1)
#> [1] "pcaCoDa"
names(prot1)
#> [1] "scores" "loadings"
#> [3] "eigenvalues" "method"
#> [5] "princompOutputClr" "mult_comp"
#> [7] "seed" "init_seed"
#> [9] "samples" "sub2D"
#> [11] "sub3D" "transformation_method"
#> [13] "transformed_data" "dist_matrix"
#> [15] "name"Simplify CoDa names
We may want to simplify the names of the transformed variables before plotting them in a biplot. The transform_coda function, which is called inside apply_ordination for protocol 1, generates composite names with format “transformationMethod-component” for all transformed variables (e.g., “CLR-Fe2O3”). The simplify_coda_names function replaces these names back to the shorter version (e.g., “Fe2O3”). However, you must always remember that the variables projected in biplots are not the originals but the transformed versions. This is particularly important when dealing with log-ratio variables since they contain information that goes beyond the original variable (i.e., divider).
prot1 <- simplify_coda_names(prot1)Test the given chemical reference groups
Perform four tests (anosim, betadisper, permdisp2, and permanova) that assess the separation and uniformity of the given group factor. For more details on these tests, we refer to:
Anderson, M.J., Walsh, D.C.I., 2013. PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions: What null hypothesis are you testing? Ecol. Monogr. 83, 557-574. doi:10.1890/12-2010.1
The whole test batch may take several minutes depending on the size of the data matrix and the number of groups.
prot1_tests <- test_groups(prot1$dist_matrix, factor_list$ChemGroup)
#> [1] "initiating test batch..."
#> [1] "vegan::anosim done."
#> [1] "vegan::betadisper done."
#> [1] "vegan::permutest done."
#> [1] "vegan::adonis done."
#> [1] "Test batch completed."The tests outputs can be accessed by their name:
names(prot1_tests)
#> [1] "permanova" "betadisp" "permdisp2" "anosim" "text"The object also contains a “text” object, which is a function that generates a list of text lines for plotting the results of PERMANOVA and PERMDISP2 tests. It may feel confusing, but keep in mind that this “portable” function requires the same ordination object as an argument.
displayTestText <- function(test_text) {
par(mar = c(0, 0, 0, 0), fig = c(0.05, 0.9, 0.05, 0.9))
plot.new()
for (i in 1:length(test_text)) {
# this is for calculating the vertical
# position of paragraphs and lines
test_spacing_paragraph = 0.8
test_spacing_line = 0.8
first_line_pos_y <-
1 - test_spacing_paragraph * ( (i - .9) / length(test_text) )
pos_y <- first_line_pos_y
if (length(test_text[[i]]) > 1) {
next_paragraph_pos_y <-
1 - test_spacing_paragraph * ( i / length(test_text) )
for (j in 2:length((test_text[[i]])))
{
pos_y <-
c(pos_y,
first_line_pos_y - test_spacing_line *
((j - 1) / (length((test_text[[i]])))) *
(first_line_pos_y - next_paragraph_pos_y)
)
}
}
# display a paragraph (a element of the list)
text(0, pos_y, labels = test_text[[i]], cex = 0.8, pos = 4)
}
}
displayTestText(prot1_tests$text(prot1_tests))
A rule-of-thumb for interpreting PERMANOVA and PERMDISP2 results is: if both p-values are low enough (e.g. < 0.05), the classification given is a good approximation of the data.
Biplots
Ordination objects are best represented in biplots, which simultaneously display the projections of observations (points) and variables (arrows) over the same space. There are several options for creating biplots in R, starting with the readily available biplot function:
biplot(prot1)
Default biplot in R
Although there are several options for customizing this default biplot function, we recommend the use of the biplot2d3d package. This package wraps a lot of possibilities in R.
library(biplot2d3d)The biplot2d3d package use functions of other packages to allow the customization of virtually all aspects of a biplot. Another important feature of this package is the creation of three-dimensional interactive biplots using the rgl package.
Biplot 2D
You can consult all tuning options available in the biplot_2d function by calling up the help file:
?biplot_2dThe default configuration will probably give you a much clearer picture than the biplot function, specially if your dataset contains more than 100 observations.
biplot_2d(prot1)
default 2D
The default setting detaches the variables projections (arrows) and places them as a miniature in the bottom-right corner. In this format they may still be interpreted, much like the North when reading a map. Remember though: the more longer arrows you see, the less each one of them is reliable when comparing point values. Here, we were ‘lucky’ for getting two nearly orthogonal variables (CaO and MgO), which means, for instance, that those observations in the top-left corner are surely more calcareous than those in the bottom-right. See the Appendix section for more details.
If detaching the arrows is not of your preference, you can disable this:
biplot_2d(prot1,
detach_arrows = FALSE,
output_type = "preview")
detach_arrows = FALSE
You can also visualize how the projection of points respond to a given typology (in this case, the chemical reference groups defined in previous studies).
biplot_2d(prot1,
groups = factor_list$ChemGroup,
output_type = "preview")
default with groups
To get a prettier plot or match your research needs, you can play with the options given by the biplot_2d.
Note that groups are by default marked using inertia ellipses. They can only be interpreted as confidence ellipses if each group can be assumed to be normally distributed in all variables considered (see more details on the scaling factor “group_ellipse_cex” in the help file). This is often reasonable concerning groups that are either too small (< 10) or contain subgroups.
In the argument “test_text” you can introduce the “text” function of the “prot1_tests” object.
biplot_2d(prot1,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0.6,
invert_coordinates = c(TRUE, TRUE),
arrow_label_cex = 0.7,
test_text = prot1_tests$text(prot1_tests),
test_cex = 0.8,
test_fig = c(0, 0.5, 0.65, .99),
output_type = "preview")
tunning appearance, groups with colors
biplot_2d(prot1,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0,
group_star_cex = 0,
group_ellipse_cex = 0,
point_type = "label",
point_label = row.names(cleanAmphorae)[!isShipwreck],
point_label_cex = 0.5,
invert_coordinates = c(TRUE, TRUE),
arrow_label_cex = 0.7,
test_text = prot1_tests$text(prot1_tests),
test_cex = 0.8,
test_fig = c(0, 0.5, 0.65, .99),
output_type = "preview")
labeled points
It is also possible to save 2D biplots into various file formats (png, tiff, jpeg, eps):
# better PNG version
biplot_2d(prot1,
ordination_method = "PCA",
invert_coordinates = c (TRUE,TRUE),
grid_cex = 2.5,
ylim = c(-.1,.1),
point_type = "point",
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 1.5,
arrow_label_cex = 2,
arrow_cex = 0.2,
arrow_lwd = 2.5,
arrow_fig = c(.6,.95,0,.35),
subtitle_cex = 2.5,
test_text = prot1_tests$text(prot1_tests),
test_fig =c(0, 0.5, 0.62, .99),
test_cex = 2,
fitAnalysis_fig = c(0,.7,.05,.5),
# saving settings
file_name = "Prot1_Biplot2D",
directory = directories$prot1,
width = 1000, height = 1000,
output_type = "png")Biplot 3D
Most 2D options are also available when generating 3D biplots. Consult the help file for details.
?biplot_3dbiplot_3d(prot1,
ordination_method = "PCA",
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
point_type = "point",
group_representation = "stars",
star_centroid_radius = 0,
star_label_cex = .8,
test_text = prot1_tests$text(prot1_tests),
test_cex = 1.25,
test_fig = c(0, 0.5, 0.65, .99))You can save an animated GIF and a PNG snapshot using the animation function.
biplot2d3d::animation(directory = directories$prot1,
file_name = "Prot1_Biplot3D")You will need to install ImageMagick to be able to generate the GIF animation.
In these images, you get what you would see when running the biplot_3d function in a regular R session.
Prot1_Biplot3D.gif
NOTE: Animated GIF will not be displayed in the pdf version of this document.
Comparing CoDa transformations
There is a lot of debate on which transformation is useful–or even valid–for analyzing geochemical compositions in Archaeometry. We show here how you can compare the results of applying different transformations to the same dataset.
First, create different ordination objects for each type of CoDa transformation that you wish to compare:
prot1_std <- apply_ordination(cleanAmphorae[!isShipwreck,],
"1", # select protocol 1
coda_override = chemVars16,
coda_transformation = "std")
#> [1] "56.05% of variance explained in 2D"
#> [1] "64.91% of variance explained in 3D"
#> [1] "Protocol 1 ended."
prot1_log <- apply_ordination(cleanAmphorae[!isShipwreck,],
"1", # select protocol 1
coda_override = chemVars16,
coda_transformation = "log")
#> [1] "62.66% of variance explained in 2D"
#> [1] "72.53% of variance explained in 3D"
#> [1] "Protocol 1 ended."
prot1_ALR <- apply_ordination(cleanAmphorae[!isShipwreck,],
"1", # select protocol 1
coda_override = chemVars16,
coda_transformation = "ALR",
# this is the divisor component
coda_alr_base = "Fe2O3")
#> [1] "70.3% of variance explained in 2D"
#> [1] "81.75% of variance explained in 3D"
#> [1] "Protocol 1 ended."
prot1_CLR <- apply_ordination(cleanAmphorae[!isShipwreck,],
"1", # select protocol 1
coda_override = chemVars16,
coda_transformation = "CLR")
#> [1] "64.16% of variance explained in 2D"
#> [1] "75.57% of variance explained in 3D"
#> [1] "Protocol 1 ended."Then, create the respective biplots:
biplot2d3d::biplot_2d(prot1_std,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0.6,
arrow_label_cex = 0.7,
ylim = c(-0.9, 0.6),
x_title = "Scaled and Centred",
x_title_cex = 1.5,
x_title_fig = c(0.05, 0.9, 0.85, 1),
output_type = "preview")
biplot2d3d::biplot_2d(prot1_log,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0.6,
arrow_label_cex = 0.7,
ylim = c(-0.6, 0.6),
x_title = "Log-scaled",
x_title_cex = 1.5,
x_title_fig = c(0.05, 0.9, 0.85, 1),
output_type = "preview")
biplot2d3d::biplot_2d(prot1_ALR,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0.6,
arrow_label_cex = 0.7,
ylim = c(-0.6, 0.6),
x_title = "Additive log-ratio",
x_title_cex = 1.5,
x_title_fig = c(0.05, 0.9, 0.85, 1),
output_type = "preview")
biplot2d3d::biplot_2d(prot1_CLR,
groups = factor_list$ChemGroup,
group_color = color_list$ChemGroup,
group_label_cex = 0.6,
arrow_label_cex = 0.7,
ylim = c(-0.5, 0.4),
x_title = "Centred log-ratio",
x_title_cex = 1.5,
x_title_fig = c(0.05, 0.9, 0.85, 1),
output_type = "preview")